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Abstract. We present the results of numerical simulations of the spherically symmetric gravitational collapse of 
^ ■ supermassive stars (SMS). The collapse is studied using a general relativistic hydrodynamics code. The coupled 

' system of Einstein and fluid equations is solved employing observer time coordinates, by foliating the spacetime 

, by means of outgoing null hypersurfaces. The code contains an equation of state which includes effects due to 

radiation, electrons and baryons, and detailed microphysics to account for electron-positron pairs. In addition 
energy losses by thermal neutrino emission are included. We are able to follow the collapse of SMS from the onset 
' ■ of instability up to the point of black hole formation. Several SMS with masses in the range 5 x IO^Mq - lO^M© 

are simulated. In all models an apparent horizon forms initially, enclosing the innermost 25% of the stellar mass. 
From the computed neutrino luminosities, estimates of the energy deposition by z/£'-annihilation are obtained. 
Only a small fraction of this energy is deposited near the surface of the star, where, as proposed recently by Fuller 
Cl|' & Shi (1998), it could cause the ultrarelativistic flow believed to be responsible for 7-ray bursts. Our simulations 

^ , show that for collapsing SMS with masses larger than 5 x 10^ Mq the energy deposition is at least two orders 

5_j ' of magnitude too small to explain the energetics of observed long-duration bursts at cosmological redshifts. In 

I addition, in the absence of rotational effects the energy is deposited in a region containing most of the stellar 

. mass. Therefore relativistic ejection of matter is impossible. 
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1. Introduction have essentially been proposed, relying upon either stel- 
lar dynamics in a dense star cluster or on the hydrody- 

Recent observations with the Space Telescope Imaging ^^^.^^ ^ ^.^^^^ supermassive star (SMS). Both cases 

Spectrograph on board the HST have added stronger sup- ^^^^ ^^^^^^^ ^ dynamical instability (Chandrasekhar 1964; 

port to the mcreasing behef that supermassive black holes p^^j^^ ^gg^. ^^^^.^^ ^ Teukolsky 1985) and can therefore 

(SMBH) are not unusual in nature: the central objects 1 -1. j- i,- t ^-l v.- v. u 

; ' 1-11 • 1 1 • undergo catastrophic gravitational collapse which would 

in more than 45 galaxies have been examined and in 34 , , , , , r ,• r oiv/mTT 

° lead to the lormation ot a SMBH. 
of them the presence of a SMBH has been confirmed 

(Kormendy 2000; see also Rees 1998 for an overview of 

SMBH). The inferred masses of these objects imply that Respite of missing theoretical and observational evi- 

SMBH of 106 - 10^ Mq are present in the center of proba- 'ience for the existence of SMS, it is still of theoretical 

bly most, if not all, galaxies. Our own Galaxy with a cen- interest to study their properties, due to their potential 

tral mass between 2.6 x IO^Mq and 3.3 x IO^Mq (Genzel astrophysical relevance as progenitors of SMBH. In par- 

et al. 2000) lies at the lower limit of this mass range. In the ^i^ular, it is worth analyzing the dynamics involved in 

particular cases of NGC 4258 and our Galaxy, the radius ^^e gravitational collapse of such stars. The low frequency 

inferred for the location of the "dark mass" is sufficiently gravitational waves emitted during such (in case of rotat- 

small to exclude any black hole alternative (Maoz 1998). st^^^s) events could be detected by the proposed Laser 

The question of how SMBH form, is, however, still not Interferometer Space Antenna (LISA). Strong bursts of 

settled and the nature of their progenitors is still rather gravitational radiation would only occur for homologous 

uncertain (see, e.g., Rees 1984). Two different scenarios collapse down to small radii, but not in case of the forma- 

tion of an initially small black hole which gradually grows 

Send offprint requests to: Jose A. Font through subsequent accretion (Thorne & Braginsky 1976). 
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A number of collapse simulations of SMS have been 
attempted through the years. The first ones were done by 
Appcnzcllcr & Fricke (1972a, b) who studied the collapse 
of spherical stars of masses 7.5 x lO^M© and 5.2 x lO^M©, 
using a code with only a subset of the general relativistic 
equations implemented. Since core; temperatures higher 
than 5 x 10^ K were encountered, nuclear burning was ex- 
pected to release high amounts of energy. They also found 
that for stars with masses greater than IO^Mq thermonu- 
clear reactions have no major effect on the evolution. In a 
subsequent paper Fricke (1973) studied the fate of SMS de- 
pending on the initial metallicity and mass. Later, Shapiro 
& Teukolsky (1979) simulated the evolution of a IO^Mq 
star with a spherically-symmetric relativistic code which 
did not include any niicrophysics effects. They were able 
to follow the whole evolution of the collapsing star until 
the formation of a black hole. Fuller, Woosley & Weaver 
(1986) revisited the work of Appenzeller & Fricke per- 
forming simulations in the range of 10^ — lO^M© with a 
code which included post-Newtonian corrections and de- 
tailed microphysics. They found that SMS with zero initial 
metallicity do not explode, but the critical metallicity for 
an explosion was considerably reduced compared to the re- 
sults of Appenzeller & Fricke. In a scries of recent papers, 
Baumgarte & Shapiro (1999a,b) investigated the influence 
of rotation on the evolution and luminosity of a SMS, find- 
ing that the luminosity is reduced by rotation and, hence, 
the lifetime of the star is enhanced. Furthermore, they 
showed that the ratios R/M, T/\W\ and J/M^, where 
R is the radius, M the mass, T the rotational energy, J 
the angular momentum and W the gravitational potential 
energy of the star, are universal numbers at the onset of 
instability. Therefore they concluded that the collapse of 
SMS should produce universal gravitational waveforms. 

SMS have recently re-gained theoretical interest since 
Fuller & Shi (1998) proposed a scenario to explain the ori- 
gin of cosmological 7-ray bursts (GRBs) from the collapse 
of SMS. The large amounts of emitted thermal neutrinos 
and antineutrinos during the catastrophic collapse of such 
stars (or, alternatively, of stellar clusters) could deposit 
a significant fraction of their energy by annihilation to 
electron-positron pairs near the surface of the star. They 
argued that this annihilation-induced heating could lead 
to relativistic expansion and associated 7-ray emission by 
cyclotron radiation and/or by the inverse Compton pro- 
cess, producing the required energies for a GRB in the 
range 10^^ — lO^** ergs within several seconds (assuming 
isotropic emission). 

In order to give reliable, quantitative estimates of the 

neutrino emission from the collapse of SMS we have per- 
formed a number of numerical simulations of spherically 
symmetric collapsing SMS, using a relativistic hydrody- 
namics code (Papadopoulos & Font 2000a). The code is 
able to follow the collapse from the onset of instability 
up to the point of black hole formation. This was an im- 
portant consideration since, due to the large temperature 
dependence of the neutrino emission, most of the neutri- 



nos are emitted close to the onset of black hole formation, 
when the temperatures in the stellar core are highest. 

This paper is organized as follows: in Section 2 we 
present the problem setup, describing the equation of state 
we use, the equations of general relativistic hydrodynamics 
and the initial data for SMS. We also discuss some com- 
putational issues relevant to our simulations. These simu- 
lations are presented and analyzed in Section 3, together 
with some tests to show the accuracy of the numerical 
code. Moreover, this section contains our results concern- 
ing the neutrino emission and the energy deposition by 
neutrino-antineutrino annihilation during SMS collapse. 
In Section 5 we discuss the implications of our results for 
GRBs. The paper ends with a summary in Section 6. 

2. Problem setup 

2.1. Properties of SMS and microphysics 

SMS are supported against gravitational collapse by ra- 
diation pressure, the gas pressure being negligible. Such 
stars are very well approximated in equilibrium by poly- 
tropes, p = Kp^, r = 1 -|- 1/n with polytropic index n = 3 
(F = 4/3). 

SMS have a two-dimensional parameter space: the cen- 
tral density pc, which controls the radius of the star, and 
the mass of the star M, which determines the ratio (3 
between gas pressure and total pressure via Eddington's 
quartic equation (Kippenhahn & Wigert 1994) 
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where p. is the mean molecular weight. The polytropic 

constant K is related to (3 by 
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where TZ is the universal gas constant and a is the radia- 
tion constant. 

The adiabatic index is slightly larger than 4/3 due to 
the contribution from baryons. The onset of instability 
occurs when the adiabatic index drops below a critical 
value 

4 ^ ,„2GM 
3+^-^2^' 
where G is the gravitational constant and c is the speed 
of light. This corresponds to a critical central density 
(Chandrasekhar 1964) 
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Stars with a central density pc larger than pcrit will become 
dynamically unstable. 

In our code we have implemented a tabulated equation 
of state (EOS) p = p{p, e) which includes radiation, elec- 
trons and effects associated with the creation of electron- 
positron pairs. Contributions due to the presence of an 
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Fig. 1. In the 3+1 formulation initial data are constructed 
on spacelike hypersurfaces, t = const. In the character- 
istic formulation employed in this work, initial data are 
constructed on null hypersurfaces (lightcones), u = const. 



additional Boltzmann gas (H+) are taken into account in 
an appropriate way. Given e, the specific internal energy, 
and /9, it is possible to determine the temperature T by a 
Newton-Raphson iteration and then compute the pressure 
from p and T . 

Furthermore, the EOS allows us to compute neutrino 
energy loss rates due to pair annihilation, + e+ ^ + 
F, photo-neutrino emission, 7 -I- — > + + F, and 
plasmon decay 7 ^ + 17. The fitting formulas and tables 
for the neutrino energy loss rates of the first two processes 
have been taken from Itoh et al. (1996). The corresponding 
expressions for plasmon decay were obtained from Haft et 
al. (1994). These processes become relevant in the very 
last epochs of the collapse. 

2.2. General relativistic hydrodynamics 

The formulation and implementation of the Einstein and 
hydrodynamic equations follows the work of Papadopoulos 
and Font (2000a, b). These equations are formulated 
adopting an outgoing null foliation of the spacetime 
(Bondi et al. 1962; Sachs 1962), whose line element, in 
spherical symmetry reads: 

""^"^ ' " - - r2(d02^sin2M02). (5) 



The radial coordinate r is chosen to make the spheres 
of rotational symmetry have an area of Anr^. The time 
coordinate represented by u will be referred to as the 
"observer time coordinate" . This retarded time u labels 
the null cones in terms of proper time along the central 
geodesic. Additionally, 6 and (j) s-re angular coordinates 
for the null rays, as shown in Figure I. The geometry is 
completely described by the two metric functions V{u, r) 
and f3{u, r). 

The hydrodynamic equations are cast into a first-order, 
flux-conservative system which, in spherical symmetry, 
reads: 

duV + -Ls.V^F = S, (6) 



Table 1. Parameters of the 
initial models: Mass and 
central density. The crite- 
rion for the onset of dy- 
namical instability, Eq. (|^), 
gives the central densities. 
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where T'"' is the stress-energy tensor, are the 

Christoffel symbols and ^/—g is the volume element, which 
is given, for the chosen geometry, by \/^g — r^e^^ sir? 6. 
The state vector - the conserved quantities - contains the 
relativistic densities of mass, momentum and energy, de- 
fined, in terms of the primitive variables, w = (p, u'",e), 
as £> = S"" = phu°u'' + pg°'' and E = phu°u° +pg°°, 
respectively. In this expressions h is the specific enthalpy, 
h = 1 + e + p/ p, and and are the time and radial 
components of the 4-velocity, respectively. 

The Einstein equations reduce to two radial hypersur- 
face equations for f3 and V: 



(10) 
(11) 



The integration of the metric is thus completely separated 
from the recovery of the primitive quantities, which sim- 
plifies the recovery. Instead of a three-dimensional root 
finding procedure, coupled to the two ODEs for the met- 
ric, only a one dimensional problem must be solved. 

2.3. Initial data 

All SMS models in our sample are chosen such that the 
initial central density is a bit larger than the critical cen- 
tral density given by Eq. (^), for the particular mass of 
the star. They are listed in Table |l|. The initial models for 
the simulations are therefore described by one parameter 
only, the mass of the star. 

The equations describing the hydrostatic equilib- 
rium are the Tolman-Oppenheimer-Volkoff (TOV) equa- 
tions. Using outgoing Bondi-Sachs coordinates they read 
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(Papadopoulos & Font 2000a): 

P,r =(^^-^{l + 87rr'p?jph, (12) 

= l + 87rr^{p- ph), (13) 

where Y = Ve~^^. These coupled differential equations 
are integrated by a fourth-order Runge-Kutta scheme. The 
relation between pressure and density is prescribed by the 
polytropic structure relation. The integration of Eq. (|lO|) 
yields /3. The boundary conditions at r = are p ~ Pc — 
Kpt^^ ,Y — and /3 = 0. The polytropic constant can be 
computed using Equation (^ and the corresponding (3 is 
given by Equation (|l|). 

It may be criticized that initial data for an unstable 
star be prescribed on a lightcone rather than on a spacehke 
hypersurface. For stable stars the density profile is inde- 
pendent of the spacetime slicing but in dynamic situations 
this does not hold anymore. Other authors (Baumgarte 
et al. 1995) circumvent this problem. They generate ini- 
tial data on null hypersurfaces by using a spacelike 3+1 
code. They trace outgoing lightrays, and store the hydro- 
dynamic quantities encountered along its path. However, 
as the stellar structure does not change much during the 
first light-crossing time, a direct prescription on null hy- 
persurfaces seems to be sufficient and justified. 

2.4. Computational issues 

The numerical code uses a conservative Godunov-type 
scheme with an approximate Riemann solver to integrate 
the (hyperbolic) hydrodynamic equations. The hypersur- 
face equations (Eqs. (|l^) and ([ll|)) for the metric compo- 
nents are solved using a two-step Runge-Kutta scheme. 
Specific details about the numerical algorithms can be 
found in Papadopoulos & Font (2000a). The canonical 
Eulerian grid used for our collapse simulations of SMS, 
consists of 500 fixed zones distributed in a non-equidistant, 
geometrically growing way from the origin, such that the 
final Schwarzschild radius of the star is resolved with 
roughly 50 zones. 

It is worth discussing in some detail the computational 
difficulties arising from the use of outgoing lightlike space- 
time foliations when a black hole forms. We have found 
that in numerical simulations of collapsing SMS (see next 
Section) a black hole forms at the center which only en- 
compasses the innermost ^25% of the total mass. The 
rest of the star's mass accretes on this seed black hole 
in a dynamical timescale. When the computation is con- 
tinued beyond the moment when the black hole first ap- 
pears, the density profile still looks reasonable, but the 
internal energy develops spikes and grows rapidly in the 
region M{r) > 0.25M. As expected, in observer time co- 
ordinates the region interior to that mass does not evolve 
further and thus remains intact. The unphysical behavior 
can be postponed by severely limiting the time step near 
black hole formation. Similar observations were reported 
by Gomez et al. (1996) who simulated the evolution of a 



scalar field in observer time coordinates up to the point of 
black hole formation. 

The conditions arising at the onset of black hole forma- 
tion are shown in Figure ^ for the collapse of a pressureless 
(dust) star of homogeneous density. Lightrays, whose tra- 
jectories can be computed analytically, are emitted from 
the center of the star. They cross the surface of the col- 
lapsing star when the formation of an apparent horizon is 
about to begin and thus they suffer a severe delay. 

The mass m of the accreting mass shell is only impor- 
tant for determining the final event horizon, its trajectory 
being independent of this mass. To approach the inner- 
most apparent horizon is most critical since the "distance" 
between the apparent horizon and the computational null 
hypersurface is smallest there. Our numerical code ad- 
vances data on outgoing null hypersurfaces. In Figure ^ 
the three trajectories of the lightrays represent the last 
three time steps of the simulation. The region below the 
uppermost lightray trajectory, particularly the segments 
of the stellar surface and of the mass shell, drawn by solid 
lines, have already been computed. Without taking addi- 
tional measures it is not possible to compute the dotted 
section of the accreting mass shell. 

Other authors have reported similar problems with 
their implementation of the observer time method in the 
context of the collapse of neutron stars to black holes 
(Baumgarte et al. 1995). In their case it helped to reduce 
the accuracy of their numerical method by performing a 
first-order accurate integration of the metric component 
goQ. Effectively, this first-order scheme constantly under- 
estimates goo I thus preventing the penetration of the hori- 
zon. In our formulation of the hydrodynamic equations 
it turns out that manipulating the metric while keeping 
the conserved quantities fixed influences the distribution 
of the primitive quantities in such a way that the met- 
ric component goo is less underestimated or even overesti- 
mated. 

Our result that an apparent horizon forms inside the 
collapsing SMS at about M{r) ^ 0.25A/sms agrees quite 
well with the results of other authors. Woosley, Wilson and 
Maylc (1986) found a trapped surface at M(r) ~ 0.2Msms 
for a 5 X IO^Mq SMS and similar results were reported 
by Shapiro & Teukolsky (1979) for a IO^Mq SMS. The 
slightly smaller value of the mass inside the apparent 
horizon most likely arises from the use of different mi- 
crophysics. 

Our main aim is to compute the lightcurves of neu- 
trino emission during the collapse of a SMS. The neutrino 
luminosities will, like lightrays, suffer from severe redshift 
when black hole formation is in progress. However, as we 
discuss later, almost all neutrinos are emitted near the 
center of the star. Therefore, the dominant part of the 
redshift for these neutrinos arises from the innermost ap- 
parent horizon - the matter further outside has not yet 
collapsed far enough, i.e., the outermost apparent horizon 
is too distant to have a significant impact on the redshift. 
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Fig. 2. Spacetime diagram of the collapse of a homoge- 
neous dust sphere in arbitrary units. The vertical axis 
shows the advance of local proper time T(r), i.e., the 
time measured by observers at fixed radius r. The tra- 
jectories of three outgoing lightrays (emitted at r = 
1.3978, 1.4278, 1.4578) are shown as they escape from the 
forming black hole. Due to the gravitational distortion of 
the spacetime, lightrays are severely delayed close to the 
apparent horizon. Subsequently, a shell of mass m is ac- 
creting onto the seed black hole. In such a situation most 
evolution codes are likely to experience problems at the 
innermost apparent horizon due to finite numerical preci- 
sion. Without taking additional measures it is not possi- 
ble to compute the dotted section of the accreting mass 
shell. The dashed sections of the spacetime diagram are 
already inside the event horizon and are causally discon- 
nected from distant observers. 



3. Simulations 

We turn next to describe our simulations, presenting first 
a series of tests we have performed in order to calibrate 
the code. 



3.1. Tests of the numerical code 

For testing the accuracy of the code we use spherical neu- 
tron star configurations as there exist either previous re- 
sults or analytic estimates that we can use. Convergence 
tests show that the code is second order accurate in space 
and time except at the center of the star. As the center is 
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Fig. 3. Test case for determining the onset of instability in 
a neutron star. The upper panel shows the mass against 
the central density pc. Models with pc > 3.96 x 10^^ g 
cm^'^ should be dynamically unstable. In the lower panel 
the surface radii is plotted as a function of time for two 
models with only slightly different central density. While 
the model with pc — 3.868 x 10^^ g cm~'^ is stable and its 
radius is kept constant in time (solid line), the one with 
Pc = 3.869 X 10^^ g cm~^ collapses to a black hole in less 
than 1 ms. 



a local extremum, the numerical schemes we use (which 
belong to the so-called total variation diminishing class) 
are only first order accurate there. The code can keep a 
polytropic model of a neutron star in equilibrium for very 
long-term evolutions (roughly 200 light-crossing times). 

Determining the onset of instability of stellar mod- 
els by probing different central densities provides a very 
strong test for numerical codes. In order to do this a se- 
quence of neutron star equilibrium configurations is con- 
structed by solving the TOV equations for several central 
densities. Only those configurations where > are 
dynamically stable. In order to benchmark our results we 
use the same model as Baumgarte et al (1995), namely a 
5/3 polytrope with K = 5.38 x 10^ (in cgs units) with a 
mass of 0.79 Mq. In Figure || the dependence of M on pc 
is shown (upper panel) . The maximum mass is obtained at 
Pc — 3.96 X lO^^g cm~'^ . The lower panel in Figure || shows 
the surface radii of two different models having almost the 
same central density. The model with pc = 3.868 x lO^^g 



cm~'^ remains stable (solid line), while the model with 
Pc = 3.869 X lO^^g cm~^ collapses (dashed line). This is 
exactly the same numerical value quoted by Baumgarte 
et al (1995). Although the determination of the critical 
density is only accurate at 2.3%, the corresponding mass 
is identical up to the first three significant figures. 

The accuracy of the code has also been measured by 
computing the frequencies of the radial modes of pulsa- 
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tions of a spherical neutron star. The numerical evolution 
of an initially static star is influenced by the local trunca- 
tion error of the hydrodynamic scheme. These errors ex- 
cite low amplitude radial pulsations whose frequencies can 
be measured by Fourier transforming the time evolution 
data of any fluid variable. Since in our code these frequen- 
cies are measured in observer time coordinates (and not 
in terms of local proper time) one has to correct them by 
a redshift factor. The results can be compared to compu- 
tations of linear normal modes (as an eigenvalue problem) 
or with results of other nonlinear codes. Our results agree 
very well with the frequencies computed by Font et al. 
(2000). Using 200 equidistant radial zones, the fundamen- 
tal radial mode and the first two harmonics are computed 
with an accuracy of 0.12%, 3% and 0.1%, respectively. 

3.2. Evolution of collapsing SMS 

We now describe the dynamics of the collapse of a spher- 
ical supermassive star. Although we have computed the 
collapse of all models listed in Table |l|, we focus our dis- 
cussion on the 5 x 1O^M0 star, which represents our canon- 
ical model. Since all SMS are very closely approximated by 
4/3-polytropes, the evolution of stars of different masses 
involves different length and time scales, but it is other- 
wise very similar. 

Figure ^ shows the evolution of the density, tempera- 
ture, metric components and velocity for the 5 x 1O^M0 
SMS. The collapse lasts 8 x lO^^ s (« 9.3 days) and the 
central density increases by a factor of 1.08 x lO"^. The 
initial configuration can be well described by Newtonian 
theory, as goo{R) = —1-0058 at the surface of the star 
R, which deviates from the value in Minkowski space- 
time, goo = —1, by only 0.58%. At the end of the sim- 
ulation the final configuration has become highly rela- 
tivistic and goo{R) = —119. According to the relation 
dr^ ~ {e^^V/r)du^ , this corresponds to a time dilation 
factor of about 11. In order to visualize this dramatic in- 
crease. Figure ^ displays the deviation from flat spacetime 
logarithmically. In all simulations the center of the star is 
chosen to resemble Minkowski geometry, e.g. the bound- 
ary condition for /? is chosen as /3(r = 0) = 0. Due to the 
use of observer time coordinates all hydrodynamic quan- 
tities seem to freeze at the end of the simulation while the 
metric quantities start to grow very rapidly when black 
hole formation is about to begin. The maximum radial 
velocity encountered is about —0.57c at the very end of 
the simulation. The velocity profile is proportional to r up 
to r w 5.6 X 10^° cm, which corresponds to about 25% of 
the star's mass. 

The panel on the lower right side of Figure ^ shows 
the local proper time r(r) versus the position of mass 
shells enclosing fixed fractions of the total mass of the star. 
Between two lines the enclosed mass increases by 10% of 
the total mass of the star. The three lines intersecting the 
trajectories of the mass shells represent hypersurfaces of 
constant coordinate time u. They can be interpreted as 



the trajectories of outgoing lightrays. The arrow indicates 
the slope of a lightray at asymptotically large radii, i.e., in 
Minkowski spacetime. Thus, lightrays are severely delayed 
due to the distortion of the spacetime. By comparing the 
position of the mass shell enclosing 90% of the mass of the 
star with the Schwarzschild radius of a 5 x IQ^Mq star, 
Rs = 1.47 X 10^^ cm, it becomes clear that the star still 
has to contract to less than one sixth of its current radius. 

The density profiles in Figure ^ are self-similar, i.e., 
the collapse proceeds homologously. A perfect homologous 
collapse is found until the central density has increased 
by a factor 798 corresponding to a central density pc w 18 
g/cTo? . At pc ~ 10'^ g/cm^ the deviation from homology is 
evident. Goldrcich & Weber (1980) found analytically that 
in Newtonian theory a star with polytropic index F = 4/3 
should contract strictly homologously. The deviation from 
homologous collapse observed in the simulation is mainly 
due to effects of general relativity: time dilation results in 
a slower advance of local proper time T(r) near the cen- 
ter of the star, and the proper volume element ^ —gdr is 
larger than in Newtonian gravity. Therefore the density 
profiles in Figure ^ tend to become steeper. Another con- 
tribution comes from the use of detailed microphysics. Due 
to copious e^e+ pair production at the center of the star, 
where temperatures are highest, the increase of the inter- 
nal energy is used to create the rest mass of the leptons 
instead of generating additional pressure. The EOS is thus 
softened, i.e. the adiabatic index is reduced, and the star 
tends to become more centrally condensed. However, since 
the deviation from homology is already apparent when the 
central temperatures are lower than lO^K (see Figure ||), 
e~e+ pair creation can be considered as a minor contri- 
bution. 

The evolution of SMS of different masses is found to be 
qualitatively very similar. The mass essentially affects the 
evolution only as a scaling parameter, except for the effects 
of neutrino emission. In all models of our sample, a black 
hole forms from the innermost 25% of the total stellar 
mass, and thus, the location of the innermost apparent 
horizon is proportional to M. In addition, the velocity 
profiles at the onset of black hole formation are found to 
be nearly identical. The timescale Ar of neutrino emission, 
Ar « Rs/u^{Rs), is thus proportional to the stellar mass 
M (see also Figure |^). The dependence of the final core 
temperature on the mass is found to be Tc ~ M""-^ (Shi 
& Fuller 1998). 

3.3. Neutrino emission 

As already mentioned, three different kinds of thermal 
processes of neutrino production are included in the simu- 
lations: e^e~-pair annihilation, photo-neutrino emission, 
and plasmon decay. The neutrino generation rates in ther- 
modynamic equilibrium depend only on the temperature 
and density of the matter. Figure |5| shows the evolution- 
ary tracks, including effects due to electron-positron pair 
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Fig. 5. Evolutionary tracks of several simulated SMS in 
the p — T-plane. More massive stars develop lower final 
core temperatures at the onset of black hole formation. 



creation in the stellar gas at T > lO^K, for our sample of 
stars in the range of 5 x lO^M© - IO^Mq. 

From Figure ^ it is obvious that the plasmon decay 
process is negligible for the simulations. For stars with 
masses below 5 x 10* Mq the pair creation process will 
be the most important process, as most of the energy 
is released in the last decade of central density increase. 
Black hole formation prevents stars with masses larger 
than 10^ Mq from entering the region dominated by pair 
annihilation. Neutrinos are then emitted primarily by the 
photo-neutrino process. Since the neutrino emission rate 
per volume for pair annihilation is Q^^ ~ [erg/s/cm"^] 
(Itoh et al. 1996), one can already infer from Figure || that 
the total energy release in form of neutrinos decreases with 
increasing stellar mass. 

As the volume increases r^dr) and the tempera- 
ture decreases with radius, most neutrinos are emitted 
from a narrow spherical shell deep inside the star. Figure 
^ depicts the radial dependence of the differential neu- 
trino luminosities, dL,y-^{r)/dr = Anr^Q^, for four models 
of our sample at various epochs (represented by the differ- 
ent lines for each star) . The shells are centered at a radius 
at which the differential luminosity has a maximum. 
The total (unredshifted) neutrino luminosity for all flavors 



Lt7., + LlJ 



is given by 



Aiir'^Q^ (T(r,t))dr. 



(14) 



(15) 



Figure |7| shows the computed neutrino luminosities 
LijT7{too) as seen by an observer at infinity for the collapse 
of different SMS. For each stellar mass the bare and the 
redshifted neutrino luminosities are depicted. It is pos- 
sible to identify the moment in time when gravitational 
redshift overcompensates the increasing neutrino energy 
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Fig. 6. Radial dependence of the (scaled) differential neu- 
trino luminosities dL,^-i^{r)/dr for four SMS: (A) 5 x 
10^ A/©, (B) 10^ Mq, (C) IO^Mq and (D) IO^Mq. The 
various lines indicate different epochs for each model (in- 
creasing upwards). The upper curves correspond to the 
end of the simulations. The scale factors used are 5 x 10''^ 
(A), 10'''' (B), 2 X (C) and 2 x 10^6 (D), respectively. 
The main fraction of the neutrinos is emitted around a 
radius R^, which is much smaller than the radius of the 
star at the end of the simulation (~ 1.47 x lO^^cm for the 
5 X IO^Mq SMS). 



release rates caused by growing temperatures. At the end 
of the simulation the luminosities are redshifted by more 
than two orders of magnitude. 

The total energy release in form of neutrinos during 
the collapse is shown in Figure ^ as a function of the stel- 
lar mass. The data points can be approximated by two 
different power laws. Doppler shift and gravitational red- 
shift reduce the plotted values by approximately a factor 
of five. We find that for a 5 x IO^Mq SMS the total radi- 
ated energy amounts to 3.0 x 10^^ ergs, corresponding to 
an efficiency for converting rest mass energy to neutrino 
energy of 3.4 x lO"**. This efficiency is depicted in Figure ^ 
and defined as the (unredshifted) energy released in neu- 
trinos divided by the total rest-mass energy of the star. 
After the formation of the black hole and after we stopped 
our simulations, we do not expect a significant fraction of 
the neutrino energy to be emitted. The reasons for this ex- 
pectation will be discussed in some detail below. Therefore 
the efficiency is normalized to the rest-mass energy of the 
whole star, Mc^. 

These results are in good agreement with the simu- 
lation of Woosley, Wilson & Mayle (1986), who included 
neutrino transport during the collapse of a 5 x lO^M© 
SMS with zero metallicity. They found a total energy out- 
put of 2.6 X 10^^ ergs in form of electron antineutrinos. 
As the neutrino energy release in the pair-dominated re- 
gion depends on the ninth power of the temperature, most 
neutrinos are emitted, in the case of a 5 x IO^Mq SMS, 
from the very center of the star within the last Atqo ~ 10 
seconds of the collapse before a black hole forms. Even 
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Fig. 7. Time evolution of the neutrino luminosities, bare 
(solid) and redshifted (dashed), for several stars in the 
range 5 x 10^ Mq — IO^Mq. Quantity Too represents the 
proper time for an observer at infinity. This time is 
measured relative to the overall collapse timescale tq = 
{8.05 X 105s,1.7x 10«*s, 8.05 x 10^s,3.2 x lO^s} and scaled by 
Ms = M/(10'^Mq). As expected, the maximum luminos- 
ity decreases with the mass of the star. The change of the 
slope of the neutrino light curves at luminosities ~ 10''^ 
erg/s, visible in the models with 10 ''M© and 10^ Mq, is 
due to the transition from photo-neutrino dominated to 
pair-dominated emission. The rates of energy release by 
neutrino production is much more sensitive to the tem- 
perature in the latter case. 



slightly different core temperatures at the onset of black 
hole formation can result in a large difference of the total 
energy release. 

The broken power law behavior can be understood in 
terms of simple considerations. Assuming adiabatic con- 
traction, Ti, ~ M~^^^, a typical timescale of neutrino emis- 
sion At ~ M (see Figure ^, a radiating volume with ra- 
dius Ru ^ Rs ^ M (see Figure ||) and a neutrino energy 
emissivity proportional to the ninth power of the tem- 
perature (pair annihilation) at R^, (Itoh et al. 1996), the 
dependence of the total energy release on the mass of the 
star is given by ~ M'^-^ (Fuller & Shi 1998). This de- 
pendence describes very well the results of the simulations 
in the low mass range (high core temperatures > 3x10^ 
K). Stars with masses above 5 x lO^M© reach core temper- 
atures just below 10^ K where the neutrino emissivity is 
proportional to the 20th power of the temperature because 
of the onset of e"'"e~-pair creation (Itoh et al. 1996). The 
total energy release is thus proportional to E^, ~ M~^, 
which is again in good agreement with the simulations. In 
summary, one finds for the total neutrino luminosities the 
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Fig. 8. The total energy release in form of neutrinos 
against the mass of the SMS. The symbols indicate the 
computed models. The data points can be described by a 
broken power law. 
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Fig. 9. Efficiency for converting the rest mass of a SMS 
to neutrinos. As less massive stars evolve to higher core 
temperatures, both the total energy emitted in form of 
neutrinos and the conversion efficiency, rise with decreas- 
ing mass. 



limiting relations 



nQ5 < M < 5 ^ ^q6^ 

(5 X 10^ < ^ < IQS) , 
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and for the total energy release (not redshifted) 
C (10^ < li^ < 5 X 10^) , 



M- 



(5xl0^<^<108). 



(17) 



3.4. vv-annihilation and energy deposition 

In 1998 Fuller & Shi proposed the collapse of a SMS as 
a possible mechanism to trigger GRBs. In order to check 
this idea quantitatively, we calculate in this section the 
amount of energy that is deposited by z/F-annihilation in 
our collapse models of SMS. In the following, the results 
obtained for the neutrino emission in the previous section 
will also be used to estimate the spatial distribution of the 
deposited energy. 

Our hydrodynamic simulations indicate that, provided 
nuclear burning is unimportant, the gravitational collapse 
of a SMS proceeds rather smoothly, involving no shocks. 
The collapse leads to the formation of a black hole consist- 
ing of the innermost 25% of the star's mass, onto which 
the rest of the stellar matter accretes. At the same time 
an energy up to 0(10^^^) ergs is released in form of neu- 
trinos. A small fraction of this energy is deposited within 
the star. In order to accelerate stellar matter to ultrarel- 
ativistic speeds, the energy deposition must take place in 
an almost baryon-free region, i.e. very close to the surface 
of the star, or in a baryon-poor funnel along the polar 
axis in case of a rotating SMS. The kinetic energy of the 
ultrarelativistic flow could then be converted to 7-rays by 
cyclotron radiation and/or by the inverse Compton pro- 
cess. 

Such ultrarelativistic flows from SMS might be pow- 
ered by the energy deposition due to j/I7-annihilation 
(Fuller & Shi 1998), 



(18) 



since this process is capable to deposit a large amount of 
energy even in vacuum. The duration of the strong neu- 
trino emission from the collapse of a 5 x IO^Mq SMS is 
about 30 seconds for an observer at infinity (see Figure ^ 
and thus in agreement with the duration of observed long 
GRBs (Piran 1999). 

For a neutrino-radiating sphere of radius the 
volume-integrated energy deposition rate due to uv- 
annihilation above is (Goodman et al. 1987; 

Cooperstein et al. 1987) 



fe + Ff). (19) 



with A = Ke,(^,r)Gp/9TTc. Here 

,2 /„^„2 ; 



2.066 X 10-32 

cm^/erg^ is the Fermi constant, L^^^^ and L^;^^^ are 
the neutrino and antineutrino luminosities at R^, (e^) and 
(e^) energy moments of the neutrino phase space distri- 
bution, evaluated by using a fitting formula for the pair 
neutrino spectrum as given by Shi & Fuller (1998). -ft'e.(^t,T) 



is a dimensionless constant, which, according to the stan- 
dard model of weak interactions, is given by 



67r V, 



1 



4sin^6l 



+ 8 sin 9w) for v^Ve 
-I- 8 sin'* 9w) for Vny. 



with Ow being the Weinberg angle [suiOw = 0.23). 
Assuming that most neutrinos are produced by e~e"'"-pair 
annihilation it is possible to compute the total energy de- 
position rate -B^cp'^^ 



K, 



L,. 



41.225 
5.769 ~R, 



(20) 



where is the temperature at Ri, and L^^j is the total 
neutrino luminosity. 

For the sake of simplicity and comparison with Fuller & 
Shi (1998), all effects due to relativistic redshift or general 
relativistic ray bending were neglected in the evaluation of 
the i/F-annihilation. Towards the end of our simulations, 
however, they play an increasingly important role. For the 
5 X 1O^M0 star, the main neutrino emitting region is lo- 
cated near 3x 10*°cm. According to Figure^, the gradient 
of the redshift factor is steepest around that radius, close 
to which the apparent horizon forms. This gradient will 
become even steeper and grow even faster as the collapse 
proceeds. Gravitational redshift will therefore quench the 
neutrino luminosity, also in the layer where most of the 
neutrinos annihilate. At larger radii the temperature is 
too low to allow for strong neutrino production. This will 
in particular be true during the phase when the forming 
black hole swallows the remaining '-^ 75% of the SMS, be- 
cause there is no resistance to the infalling gas, and com- 
pressional heating will therefore stay moderate. For this 
reason it is very unlikely that the accretion of the rest of 
the star's mass onto the seed black hole will yield a large 
contribution to the total neutrino emission and energy re- 
lease by ^/jz-annihilation. 

The resulting energy deposition for different SMS is 
plotted in Figure |l^. The corresponding energy deposi- 
tion efficiency E^cp/Ei, is shown in Figure ^ The depen- 
dence of the deposited energy on the mass of the star 
at both ends of the considered mass range is i?dcp ~ 
L%T^At/R^ ~ Ll-M-°-^, where At is the timescale of 
the emission. Therefore, using Eq. (WEi) 




(10^ < ^ < 5 X 10^) 
(5 X 10^ < ^ < 108) 



(21) 



Fuller & Shi (1998) gave analytic estimates for the ex- 
pected neutrino emission from the collapse of a SMS. They 
employed a simple model to deduce the final average core 
temperature Tc as a function of the mass of the homol- 
ogously collapsing core M^'~^. Assuming neutrino energy 
release rates Q,^ ~ and a radiating sphere with a radius 
equal to the Schwarzschild radius, Rg = 2GM^^/c^, for 
the homologously collapsing mass, their estimate for the 
characteristic neutrino luminosity is 



i^T7 w 5 X 10"(Af5"^)-3/2 gj.g/s^ 



(22) 
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Table 2. Comparison of the analytic model of 
Fuller & Shi (1998) with the results of our numerical 
simulation. Values are missing where the approxi- 
mations by Fuller & Shi are not applicable or do not 
provide information. Overestimating the tempera- 
ture Tn at the neutrino-radiating sphere Ri, leads 
to a large error in the rate of energy release by neu- 
trino emission Qi, and thus in the energy deposition 
i?dcp at 3i?i/ < r < oo. 
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Fig. 10. The total neutrino energy deposition as a func- 
tion of the mass of the SMS. Due to decreasing neutrino 
luminosities the total energy deposited within the star de- 
creases rapidly with the mass of the star. The dependences 
on both extremes of the considered mass range can be un- 
derstood by simple considerations (see text for details). 
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Fig. 11. Efficiency for conversion of neutrino energy to 
pair plasma by z/F-annihilation. 



where Mf^ = A/^^/IO^Mq. Furthermore, assuming a 



typical emission timescale of Ar 



for the energy deposited above a radius r is 
i?dcp ~ 2.5 X lQ''\Mf'')-^-\Rjrf erg. 



their result 



(23) 



The simulation presented in the previous section shows 
that in a 5 X 10^ Mq SMS an inner core with 25% of the 
stellar mass is collapsing homologously, i.e. M^'^ = 1.25. 



M = 5 X 10'' Mq 


Simulation 


Fuller & Shi 


r. [K] 


8.0 X 10" 


11.6 X lO** 


Rv [cm] 


3.0 X 10^° 


3.8 X 10^° 


Q„{R^) [erg/s/cm^] 


4.1 X 10^^ 


1.2 X 10^^ 


Ly-n [erg/s] 


1.7 X lO'^'^ 


3.6 X lO'^^ 


-Bdcp(T- > Rv) [erg] 


4.9 X lO'^^ 




-Edop(r > 3i?^) [erg] 


8.0 X 10''® 


4.7 X lO'^^ 


M{r > 3R,) 


~ 0.65MsMS 





Table g shows the estimates of Fuller & Shi in compar- 
ison to the results of our simulation. Since the physics 
involved in the description of SMS is well known, all their 
basic assumptions have been confirmed by our simulation: 
(1) The core temperature Tc - M'^-^, (2) the ra- 

dius of the radiating sphere Ri, ~ M, (3) the time scale 
of the emission At ^ M, and (4) the computed neutrino 
energy generation rates do not differ by more than 25^0 if 
they are evaluated for the same temperature T « 10^° K. 
The large discrepancy of the total deposited energy E^cp is 
caused by the strong dependence of the energy deposition 
rate_Edep on and R^, Edep ~ T^^Rl (see Eq. |2^ with 
Eq. 



16) 



Misjudging slightly both the radius i?^ of the 
radiating sphere and the temperature at that radius 
leads to a very large overestimation of the energy deposi- 
tion rate. Furthermore, the way in which we compute the 
neutrino luminosities in our simulations, Eq. (|l^), is also 
more accurate than the formula employed by Fuller & Shi, 
L^j; = Qv (T^) ^RI- This, however, decreases the ratio of 
Fuller and Shi's neutrino luminosity to our luminosity by 
a factor of three, corresponding to a reduction of the dis- 
agreement of the energy deposition by z/i/-annihilation by 
a factor of nine, to yield a remaining discrepancy of nearly 
a factor of 600. 



4. Implication for gamma ray bursts 



4.9x10^1 



The integral amount of deposited energy, E'dcp 
erg in the case of a 5 x IO^Mq SMS, is at the lower 
end of the observed energy released in cosmic GRBs. All 
of this energy would have to be collimated into narrow 
jets with a coUimation of Sn/An < 1/1000—1/100 in or- 
der to explain an equivalent isotropic energy of lO^'^ to 
^ 10^^ erg as inferred for observed long-duration bursts 
with known cosmological redshifts (e.g., Frail et al. 2001). 
However, 99.837% (Goodman et al. 1987) of this energy 
is deposited in a spherical layer deep inside the star at 
Ru <r < 3Ri, , and only a tiny fraction near the surface of 
the star, where excessive baryon loading could be avoided. 
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Therefore ultrarelativistic ejection of matter with Lorentz 
factors 7^1 cannot be expected in spherical models: 



7 « 1 + 



^dcp(r- > Ru) 
M{r > R^)c^ 



1. 



(24) 



We therefore conclude that energy deposition by lyV- 
annihilation in the spherical collapse of a 5 x lO^M© SMS 
does not meet the demands for being a successful central 
engine for a GRB. 

A rotating SMS would have a reduced baryon density 
along the polar axis and the time scale of neutrino emis- 
sion could be longer, because a disk stabilized by cen- 
trifugal forces would prolong the accretion of matter into 
a forming black hole. Thus, the conditions to fulfill the 
constraints of GRB models might appear more suitable 
in this respect. Nevertheless, as the energy release rates 
by neutrino emission and the energy deposition in non- 
rotating models are several orders of magnitude below the 
estimates of Fuller & Shi, it still seems unlikely that the 
total energy release in a baryon-poor environment would 
be sufficiently large in rotating scenarios. However, an ul- 
timate statement requires detailed numerical simulations 
of the gravitational collapse of rotating SMS. 

As shown in Figure |lo| the total energy deposition de- 
pends strongly on the mass -Edep ~ M~^-^. Therefore, 
more massive stars are less promising to meet the energy 
requirements of a GRB. Extrapolating the total deposited 
energy of a 5 x IQ^Mq star to less massive stars, the 
expected value for a IO^Mq SMS is E^dep « 1-3 x lO^^^ 
erg. The energy deposited within [3i?^,cx)] is a factor 
1.63 X 10"'^ smaller and therefore of the order of 10^^ 
ergs, which is at the lower end of the observed energy 
range. However the baryon loading problem is still crit- 
ical. We have not attempted to simulate the evolution 
of such a star because in this case nuclear burning be- 
comes relevant. The energy release in the form of neutri- 
nos will also depend strongly on the fate of the star: If 
the energy release due to nuclear burning cannot inhibit 
black hole formation the neutrino emission increases with 
smaller stellar mass. If nuclear energy is liberated rapidly 
enough, the star will be destroyed in a thermonuclear ex- 
plosion. Consequently, there will be much less neutrino 
emission, since the core temperature will be considerably 
lower. Only detailed numerical simulations of the collapse 
of such stars can determine the expected energy release 
and their final fate. 

5. Summary 

We have presented results of numerical simulations of the 
spherically symmetric gravitational collapse of supermas- 
sive stars. Such simulations were performed using a gen- 
eral relativistic hydrodynamics code. The coupled sys- 
tem of Einstein and fluid equations was solved adopt- 
ing a spacetime foliation with outgoing null (character- 
istic) hypersurfaces. The code includes a tabulated equa- 
tion of state which accounts for contributions from radia- 
tion, electrons, electron-positron pairs and baryonic gases. 



Energy losses by thermal neutrino emission were taken 
into account. 

We were able to follow the collapse of SMS from the 
onset of instability up to the point of black hole forma- 
tion. Several SMS with masses in the range of 5 x IO^AIq — 
IO^Mq were simulated, showing that an apparent horizon 
forms in all cases, enclosing the innermost 25% of the stel- 
lar mass. This is in good agreement with previous simula- 
tions by other authors. We did not attempt to follow the 
subsequent accretion of the remaining mass onto the cen- 
tral black hole because, due to finite numerical precision, 
the innermost apparent horizon is ultimately penetrated 
(see Figures ^ and U, in particular the little kink develop- 
ing in the lightcones in frame (f) of the latter). This effect 
however, does not interfere with the main aim of our work, 
namely the computation of neutrino luminosities, because 
most neutrinos are emitted near the very center of the 
star. Therefore we did not try possible remedies to fol- 
low the evolution for longer times, such as excising the 
inner regions of the domain, in which the evolution essen- 
tially freezes, and imposing appropriate boundary condi- 
tions there. 

The neutrino luminosities for several SMS models have 
been computed and it has been possible to understand the 
dependence of the total energy release in neutrinos, E,y, 
on the stellar mass M in terms of simple scaling laws. 
Based on these results we have given estimates of the en- 
ergy deposition by i/i?-annihilation in the star. Our sim- 
ulations show that for collapsing SMS with masses larger 
than 5 x lO^M© this energy deposition is more than two 
orders of magnitude smaller than the estimates of Fuller 
& Shi (1998). In addition, all of this energy is deposited 
deep inside the star and cannot drive relativistic outflows. 
Therefore, the spherical collapse of such a SMS does not 
meet the demands for being a successful GRB model. 
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Fig. 4. Snapshots of p, T, goo, goi and arc shown for a 5 x IO^Mq SMS. The coordinate time interval Am elapsed 
between adjacent profiles diminishes since the evolution accelerates. The coordinate time is exactly equal to the proper 
time r elapsed at the center of the star. The deviation of the metric from Minkowski spacetime, goi = —1, is shown in 
the middle panels. As it is not possible to display a zero deviation in a logarithmic plot, the deviation from goo = —0.99 
is shown. In these plots log(— ^Oft — 0.99) = —2 for /i G {0, 1} represents flat spacetime. The plot on the lower right 
side shows the location of mass shells (AM = 5 x 10** M©) versus local proper time r. The lines intersecting with the 
mass shells are hypersurfaccs of constant coordinate time u. They represent trajectories of outgoing lightrays. The 
arrow gives the slope of a lightray at r ^ oo. Gravity already causes a severe delay of outgoing lightrays. Notice the 
kink in the upper hypersurface very close to the center of the star indicating the formation of an apparent horizon. 



